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Abstract 

The paper studies the problem of recover- 
ing a spectrally sparse object from a small 
number of time domain samples. Specifically, 
the object of interest with ambient dimen- 
sion n is assumed to be a mixture of r com- 
plex multi-dimensional sinusoids, while the 
underlying frequencies can assume any value 
in the unit disk. Conventional compressed 
sensing paradigms suffer from the basis mis- 
match issue when imposing a discrete dictio- 
nary on the Fourier representation. To ad- 
dress this problem, we develop a novel non- 
parametric algorithm, called enhanced ma- 
trix completion (EMaC), based on structured 
matrix completion. The algorithm starts by 
arranging the data into a low-rank enhanced 
form with multi-fold Hankel structure, then 
attempts recovery via nuclear norm mini- 
mization. Under mild incoherence condi- 
tions, EMaC allows perfect recovery as soon 
as the number of samples exceeds the order 
of ©(y'log^n). We also show that, in many 
instances, accurate completion of a low-rank 
multi-fold Hankel matrix is possible when the 
number of observed entries is proportional to 
the information theoretical limits (except for 
a logarithmic gap) . The robustness of EMaC 
against bounded noise and its applicability to 
super resolution are further demonstrated by 
numerical experiments. 



1. Introduction 

A large class of practical applications features high- 
dimensional objects that can be modeled or approx- 
imated by a superposition of spikes in the spec- 
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tral (resp. time) domain, and involves estimation 
of the object from its time (resp. frequency) do- 
main samples. A partial list includes medical imag- 
ing (Lustig et al., 2007), radar systems (Potter et al., 
2010); seismic imaging (Borcea et al., 2002), mi- 
croscopy (Schermelleh et al., 2010), etc. The data ac- 
quisition devices, however, are often limited by phys- 
ical and hardware constraints, precluding sampling 
with the desired resolution. It is thus of paramount 
interest to reduce the sensing complexity while retain- 
ing recovery resolution. 

Fortunately, in many instances, it is possible to recover 
an object even when the number of samples is far be- 
low the ambient dimension, provided that the object 
has a parsimonious representation in the transform do- 
main. In particular, recent advances in compressed 
sensing (CS) (Candes et al., 2006) popularize the non- 
parametric methods based on convex surrogates. Such 
tractable methods do not require prior information on 
the model order, and are often robust against noise. 

Nevertheless, the success of CS relies on sparse repre- 
sentation or approximation of the object of interest in a 
finite discrete dictionary, while the true parameters in 
many application are actually specified in a continuous 
dictionary. For concreteness, consider an object x [t) 
that is a weighted sum of X-dimensional sinusoids at 
r distinct frequencies {fi G [0, l]'^ : 1 < i < r}. Con- 
ventional CS paradigms operate under the assump- 
tions that these frequencies lie on a pre-determined 
grid on the unit disk. However, cautions need to be 
taken when imposing a discrete dictionary on continu- 
ous frequencies, since nature never poses the frequen- 
cies on the pre-determined grid, no matter how fine the 
grid is (Chi et al., 2011; Duarte & Baraniuk, 2012). 
This issue, known as basis mismatch between the true 
frequencies and the discretized grid, results in loss of 
sparsity due to spectral leakage along the Dirichlet ker- 
nel, and hence degeneration in the performance of CS 
algorithms. While one might impose finer gridding 
to mitigate this weakness, this approach often leads 
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to numerical instability and high correlation between 
dictionary elements, which significantly weakens the 
advantage of these CS approaches (Tang et al., 2012). 

In this paper, we explore the above spectral compressed 
sensing problem, which aims to recover a spectrally 
sparse object from a small set of time-domain sam- 
ples. The underlying (possibly multi- dimensional) fre- 
quencies can assume any value in the unit disk, and 
need to be recovered with infinite precision. To ad- 
dress this problem, we develop a nonparametric al- 
gorithm, called enhanced matrix completion (EMaC), 
based on structured matrix completion. Specifically, 
EMaC starts by converting the data samples into an 
enhanced matrix with (/•sT-fold) Hankel structures, and 
then solves a nuclear-norm minimization program to 
complete the enhanced matrix. We show that, under 
mild incoherence conditions, EMaC admits exact re- 
covery from 0{r log" n) random samples, where r and 
n denote respectively the spectral sparsity and the am- 
bient dimension. Additionally, we provide theoretical 
guarantee for the low-rank Hankel matrix completion 
problem, which is of great importance in control, nat- 
ural language processing, computer vision, etc. To the 
best of our knowledge, our results provide the first 
theoretical bounds that are close to the information 
theoretic limit. Furthermore, numerical experiments 
demonstrate that our algorithm is robust against noise 
and is applicable to the problem of super resolution. 

2. Connection to Prior Art 

The spectral compressed sensing problem is closely re- 
lated to harmonic retrieval, which seeks to extract 
the underlying frequencies of an object from a col- 
lection of its time-domain samples. This spans many 
signal processing applications including radar localiza- 
tion systems (Nion & Sidiropoulos, 2010), array imag- 
ing systems (Borcea et al., 2002), wireless channel 
sensing (Sayeed & Aazhang, 1999; Gedalyahui et al., 
2011); etc. In fact, if the time-domain representation 
of an object can be estimated accurately, then its un- 
derlying frequencies can be identified using harmonic 
super-resolution methods. 

Conventional approaches for these problems, such as 
ESPRIT (Roy & Kailath, 1989) and the matrix pencil 
method (Hua, 1992), are based on the eigenvalue de- 
composition of covariance matrices constructed from 
equi-spaced samples, which can accommodate infinite 
frequency precision. One weakness of these techniques 
lies in that they require prior information on the model 
order, that is, the number of underlying frequency 
spikes or, at least, an estimate of it. Besides, their 
performance largely depends on the knowledge of noise 
spectra; some of them are unstable in the presence of 



noise and outliers (Dragotti et al., 2007). 

Nonparametric algorithms based on convex optimiza- 
tion differ from the above parametric techniques in 
that the model order does not need to be specified 
a priori. Recently, Candes and Fernandez- Granda 
(2013) proposed a total- variation minimization algo- 
rithm to super-resolve a sparse object from frequency 
samples at the low end of its spectrum. This al- 
gorithm allows accurate super-resolution when the 
point sources are appropriately separated, and is sta- 
ble against noise (Candes & Fernandez-Granda, 2012). 
Inspired by this approach. Tang et. al. (2012) devel- 
oped an atomic norm minimization algorithm for line 
spectral estimation from ©(r logr logn) time domain 
samples. This work is limited to 1-D frequency models 
and assumes randomness in the data model. 

In contrast, our approach can accommodate multi- 
dimensional frequencies, and only assumes randomness 
in the observation basis. The algorithm is inspired by 
the recent advances of matrix completion (MC) prob- 
lem, which aims at recovering a low-rank matrix from 
partial entries. It has been shown (Candes & Recht, 
2009; Gross, 2011) that exact recovery is possible 
via nuclear norm minimization, as soon as the num- 
ber of observed entries is on the order of the infor- 
mation theoretic limit. Encouragingly, this line of 
algorithms is also robust against noise and outliers 
(Negahban & Wainwright, 2012). Nevertheless, the 
theoretical guarantees of these algorithms do not apply 
to the more structured observation models associated 
with Hankel structure. Consequently, direct applica- 
tion of existing MC results yields pessimistic bounds 
on the number of samples, which is far beyond the 
degrees of freedom underlying the sparse object. 

The rest of this paper is organized as follows. Section 3 
describes the data model and the EMaC algorithm. 
Section 4 presents the theoretical guarantee of EMaC, 
with a proof outlined in the Appendix. Numerical val- 
idation of EMaC is given in Section 5. We then discuss 
the extension to low-rank Hankel matrix completion in 
Section 6, and conclude the paper in Section 7. 

3. Model and Algorithm 

Assume that the object of interest x (t) can be modeled 
as a weighted sum of X-dimensional sinusoids at r 
distinct frequencies fi G [0, 1]^ (1 < « < r), i.e. 



x(t) = Vrf,e^"2-<*-^-) 



i=l 



(1) 



where d^'s denote the complex amplitude. For con- 
creteness, our discussion is mainly devoted to a 2- 
dimensional (2-D) frequency model. This subsumes 
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1-D line spectral estimation as a special case, and in- 
dicates how to address niulti-diniensional models. 

3.1. 2-D Frequency Model and Problem Setup 

Consider a data matrix X — (xi.i)„^,^ r.^i^ of 
size ni x 712- Suppose each entry can he expressed as 



Xkd = '^d,y^zl 



where y^ = exp (j27r/ii) and Zi ~ exp (j27r/2i) for some 
set of frequency pairs {(/ii,/2i) | 1 < « < ?"} (normal- 
ized by the Nyquist rate). We can then write X in a 
matrix form as follows: 



X = YDZ^ 



(2) 



where D :— diag[di,d2, ■ ■ ■ ,dr], and Y and Z are 
defined as 



Y := 



Z := 



1 


1 


1 


yi 


2/2 


Vr 


y'r' 


2/2 


■ yl''- 


1 


1 


1 


Zl 


Z2 


Zr 



(3) 



^"2 — 1 „n2 — 1 



,"2-1 



(4) 



Suppose that there exists a location set il of size m 
such that Xk,i is observed iff {k, I) G ft, and assume 
that f2 is sampled uniformly at random. We are inter- 
ested in recovering X from its partial observation on 
the location set il. 

Before continuing, we introduce a few notations that 
will be used throughout. The spectral norm (opera- 
tor norm), the Frobenius norm, and the nuclear norm 
(sum of singular values) of a matrix X are denoted 
by ||X||, ||X||p and ||X||^, respectively. The inner 
product between two matrices is defined by {B, C) = 
trace (B*C). We let Vn be the orthogonal projection 
onto the subspace of matrices that vanishes outside 17. 

3.2. Enhanced Matrix Completion (EMaC) 

One might naturally attempt recovery by applying the 
matrix completion algorithms (Candes & Recht, 2009) 
on the original data matrix, arguing that when r is 
small, perfect recovery of X is possible from partial 
measurements. Specifically, this corresponds to the 
following optimization program 



M||, 
subject to Vn (M) = Vn {X) , 



minimize 

JVfgC"l>'"2 



(5) 



which is the convex relaxation of the rank minimiza- 
tion problem. However, generic matrix completion al- 
gorithms (Gross. 2011) require at least the order of 
r max (711,712) samples, which far exceeds the degrees 
of freedom (which is O {r log n) due to a coupon col- 
lector's effect) in our problem. What is worse, when 
r > min (711,712) (which is possible since r can be as 
large as 711712), X is no longer low-rank. This moti- 
vates us to construct other forms (e.g. Hua (1992)) 
that better capture the harmonic structure. 

Specifically, we adopt an effective enhanced form of 
X based on 2-fold Hankel structure as follows. The 
enhanced matrix X^ with respect to X is defined as a 
fci X (711 — fci -I- 1) block Hankel matrix 



Xr 



X, 



Xo 


Xi ■ 


Xn^-ki 


Xi 


X2 ■ 


-^ni-fci+1 


^ki-i 


Xk, ■ 


Xni-1 



(6) 



where each block is a fc2 x (712 — fc2 + 1) Hankel matrix 
defined such that for all < ? < ni: 



xi,o 



Xl,l 

a;/, 2 



XJk^-l Xlk2 



•^Ln2 — k2 
X!n2 — 1 



This form allows us to derive, through algebraic ma- 
nipulation or a tensor product approach, the Vande- 
monde decomposition of each row of X as follows 



yi{0<l< m) : Xi = ZlYJDZji, 
where .Zl, .Zr and Y^ are defined as 



(7) 



Zi. := 



1 

Zl 



1 

22 



^fe2-l ^fe2-l 



1 Zl 
1 Z2 

1 Zr 



■,k2-l 



n2-k2 
'1 

112 — ^2 



yn2 — k2 



and yd := diag [yi , ?;2 , ■ ' ' ,2/r]- Plugging (7) into (6) 
yields the following: 



Xe = 






D 



ZR,YdZji, 



Y]^i-kiz^ 



' -'d 



y/{ni-ki+l){n2-k2 + l)ER 
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where Ei^ and Eji characterize the column and row 
space of Xc, respectively. The effectiveness of the en- 
hanced form relies on the shift-invariance of harmonic 
structures (any k consecutive samples lies in the same 
subspace), which promotes the use of Hankel matrices. 

One can now see that Xc is low-rank or, rank (Xc) < 
r. We then attempt recovery via the following En- 
hanced Matrix Completion (EMaC) algorithm: 



(EMaC) 



minimize 



\M, 



(8) 



subject to Vn (M) = Vn (X) , 



which minimizes the nuclear norm of the enhanced 
form over the constraint set. This convex program 
can be solved using off-the-shelf semidefinite program 
solvers in a tractable manner. 

3.3. Higher-Dimensional Frequency Model 

The EMaC method extends to higher dimensional fre- 
quency models without difficulty. For /v-dimensional 
frequency models, one can convert the original data 
to a K-io\d Hankel matrix of rank at most r. For in- 
stance, consider a 3 dimensional (3-D) model such that 

^h,h,h ~ Si=i '^iVi^ ^i^''^t ■ "'^^ enhanced form can be 
defined as a 3-fold Hankel matrix such that 



X, 



Xi, 



^l,c 
^2,e 



^A;3-l,o X 



ki ,e 



^ns — fca+l.c 
^na — l,c 



where X.^ o denotes the 2-D enhanced form of the ma- 
trix consisting of all entries a;/ j^^/j^ig obeying Z3 = i. One 
can verify that X^ is of rank at most r, and can thus 
apply EMaC on the 3-D enhanced form Xc- To sum- 
marize, for iiT-dimensional frequency models, EMaC 
minimizes the nuclear norm over all K-io\d Hankel ma- 
trices consistent with the observed entries. 

3.4. Noisy Data 

In practice, measurements are always corrupted by a 
certain amount of noise. To make our model and al- 
gorithm more practically applicable, we can replace 
our measurements by X°i through the following noisy 
model 

yii,l)en: X°=X,i+Na, 

where X°i is the observed {i,l)-th. entry, and Nu 
denotes the noise. Suppose that the noise satisfies 
WVn (-/V)||p < S, then EMaC can be modified as: 



(EMaC-Noisy) : minimize ||Mo 



(9) 



subject to \\Vn (M - X°)llp < 6. 



4. Main Results 

Encouragingly, under certain incoherence conditions, 
the simple EMaC enables accurate recovery of the true 
data matrix from a small number of noiseless time- 
domain samples, and is stable against bounded noise. 

For convenience of presentation, we denote Clc{i,l) as 
the set of locations in the enhanced matrix X^. con- 
taining copies of Xi,i, and let Wi./ := |r2o(i,OI- For 
each (i,l), we use j4.(j ;) to denote a basis matrix that 
extracts the average of all entries in il^ {i, I), i.e. 

, if (a,/3) e rJe(i,0, 



(^(.0)„,, - [f^^ 



else. 



4.1. Incoherence Measures 

In general, matrix completion from a few entries is 
hopeless unless the underlying structure is uncorre- 
lated with the observation basis. This inspires us to 
define certain incoherence measures. Let Gl and Gr 
be r X r correlation matrices such that for any i y^ j, 



(Gl 



1 1 - {ylyjP 1 - {z*z,)'' 



( * \ni— fci+l 



1 - {ztZj) 



■^2 — ^2+1 



^ -^^"J- (ni-fci + l)(l-y*y,)(n2-fc2 + l)(l-42i)' 

with the convention (Gl)jj = (GyC)^^ — 1. Note that 
Gl and Gr can be obtained by sampling the 2-D 
Dirichlet kernel, which is frequently used in Fourier 
analysis. Our incoherence measure is defined as fol- 
lows. 

Definition 1 (Incoherence) Let X^ denote the en- 
hanced matrix associated with X, and suppose the 
SVD of Xe is given by X^ = UAV* . Then X is said 
to have incoherence (p-i, fJ.2, f-a) tf i^^V o-^^ respectively 
the smallest quantities such that 

fTmin (Gl) > l//il, (Tmin (Gr) > l//ii; (10) 

2 



max 

{i,l)e[ni]x[n2 



i^e(*,or 



< 



A*2?' 



(11) 



and 



E 

aG[ni]x[n2 



\{UU*AtVV\^UJ^Aa,)f < -(^^ujb 



ni?i2 



holds for all b G [ui] x [712] 



(12) 



Some brief interpretations of the above incoherence 
conditions are in order: 

Condition (10) specifies certain incoherence among the 
locations of frequency pairs, which does not coincide 
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with and is not subsumed by the separation condi- 
tion required in (Candes & Fernandez- Granda, 2013; 
Tang et al.. 2012). The frequency pairs can be spread 
out (e.g. when their locations are generated in some 
random fashion), or minimally separated (e.g. when 
they are small perturbation of a fine grid). 

Condition (11) can be satisfied when the total en- 
ergy of each skew diagonal of UV* is proportional 
to the dimension of this skew diagonal. This is in- 
deed a weaker condition than the one introduced in 
(Candes & Recht, 2009) for matrix completion, which 
requires uniform energy distribution over all entries of 
UV* . For instance, an ideal ^2 can often be obtained 
when the complex phase of all frequencies are gener- 
ated in some random fashion. 

Condition (12) is an incoherence measure based on the 
(A'-fold) Hankel structures. For example, one can rea- 
son that a desired fi^ can be obtained if the magnitude 
of all entries of UU* AaW* is mostly even. Condi- 
tion (10) and (12) depend only on the frequency pairs. 
In fact, one can verify that UU* = Ej, (ElEi^)'^ El 
and W* = E^ {EnEl^y'^ En, which depend only on 
the locations of the frequency pairs. Condition (11), 
however, might also rely on the amplitudes D. 

Finally, the incoherence measures {ni, fi2, fJ-s) are mu- 
tually correlated, which is supplied as follows. 

Lemma 1 The incoherence measures {fXi, fi2, fJ's) of 
Xe satisfy 



M2 < f^clr, 



and 



M3 < fJ-jclr. 



(13) 



4.2. Theoretical Guarantees 



With the above incoherence measures, the main theo- 
retical guarantee is supplied in the following theorem. 

Theorem 1 Let X be a data matrix with matrix form 
(2), and CI the random location set of size m. Define 

Cs := max {^^, -, ■, — rw^ — ; — ptt I • if oil measure- 

ments are noiseless, then there exists a constant ci > 
such that under either of the following conditions: 

• Conditions (10), (11) and (12) hold and 

m > cimax(/xiCs,/i2,/^3Cs)?'log^ {nin2) ; (14) 

• Condition (10) holds and 

m > ci/^^c^r^ log (711712); (15) 

X is the unique solution of EMaC with probability ex- 
ceeding 1 — (711712)^^. 



Note that (15) is an immediate consequence of (14) 
by Lemma 1. Theorem 1 states the following: (1) 
under strong incoherence condition (i.e. given that 
(/ii, /i2, /is) are all constants), prefect recovery is pos- 
sible as soon as the number of measurements exceeds 
the order of r log (77,1712); (2) under weak incoherence 
condition (i.e. given only that /ii is a constant), per- 
fect recovery is possible from O (r^ log (711712)) sam- 
ples. Since there are at least Q{r) degrees of freedom 
in total, this establishes the near optimality of EMaC 
under strong incoherence condition. 

We would like to note that while we assume random 
observation models, the conditions imposed on the 
data model are deterministic. This is different from 
(Tang et al., 2012), where randomness are assumed for 
both the observation model and the data model. 

On the other hand, our method enables stable recovery 
even when the time-domain samples are noisy copies of 
the true data. Here, we say the recovery is stable if the 
solution of EMaC- Noisy is "close" to the ground truth. 
To this end, we establish the following theorem, which 
is a counterpart of Theorem 1 in the noisy setting. 

Theorem 2 Suppose X° is a noisy copy of X that 
satisfies WVaiX — X°)\\y < <^. Under the conditions of 
Theorem 1, the solution to EMaC-Noisy in (9) satisfies 



IX, -X, 



cllF 



< 



2^71l7l2 + 8771772 



%^/2nl^ 



with probability exceeding 1 — (711712) ^. 

Theorem 2 basically implies that the recovered en- 
hanced matrix (which contains 0{n\n'^) entries) is 
close to the true enhanced matrix at high signal-to- 
noise ratio. In particular, the average entry inaccu- 
racy is bounded above by 0( "^"^ (5). We note that in 
practice, EMaC-Noisy usually yields better estimate, 
possibly by a polynomial factor. The practical applica- 
bility will be illustrated in Section 5 through numerical 
examples. 

5. Numerical Experiments 

5.1. Phase Transition for Exact Recovery 

We examine phase transition of the EMaC algorithm 
to evaluate its practical ability. A square enhanced 
form was adopted with 77 1 = 772, which corresponds 
to the smallest Cg. For each (r,m) pair, we generated 
a spectrally sparse data matrix X by randomly gen- 
erating r frequency spikes in [0,1] x [0,1], and sam- 
pled a subset J7 of size 777 entries uniformly at ran- 
dom. The EMaC algorithm was conducted using CVX 
with SDPT3. Each trial is declared successful if the 
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M 



- 0.9 

- 0.8 

- 0.7 



100 120 140 



m: number of samples 

Figure 1. Phase transition diagrams wliere spilce locations 
are randomly generated. The results are shown for the case 
where ni — n2 — 15. 



normalized mean squared error \\X — X||f/|1-X'||f < 
10"'', where X denotes the estimate obtained through 
EMaC. The empirical success rate is calculated by av- 
eraging over 100 Monte Carlo trials. 

Fig. 1 illustrates the results of these Monte Carlo ex- 
periments when X is a 15 X 15 matrix. The empirical 
success rate is reflected by the color of each cell. It can 
be seen that the number of samples m grows approx- 
imately linearly with respect to the spectral sparsity 
r, in line with our theoretical guarantee in Theorem 1. 
This phase transition diagram validates the practical 
applicability of our algorithm in the noiseless setting. 

5.2. Stable Recovery via Singular Value 
Thresholding 

The above experiments were conducted using the ad- 
vanced semidefinite programming solver SDPT3. This 
and many other popular solvers (like SeDuMi) are 
based on interior point methods, which are typically 
inapplicable to large-scale data. In fact, SDPT3 fails 
to handle an n x n data matrix when n exceeds 19, 
which corresponds to a 100 x 100 enhanced matrix. 

One solution for large-scale data is the first-order al- 
gorithms tailored for MC problems, e.g. the singu- 
lar value thresholding (SVT) algorithm developed in 
(Cai et al., 2010). We propose a modified SVT algo- 
rithm in Algorithm 1 to exploit the Hankel structure. 

In particular, 2?rt(') in Algorithm 1 denotes singular 
value shrinkage operator. Specifically, if the SVD of X 
is given by X = UJIV* with S = diag ({ctj}), then 

VrAX):=Udmg{{{a,-Tt)^})V*. 



Algorithm 1 Singular Value Thresholding for EMaC 
Input: The observed data matrix X° on the loca- 
tion set n. 

initialize let X° denote the enhanced form of 
Pn (X°); set Mo == X° and t = 0. 



repeat 

l)Qt< 

2)Mt 
3) t 



Vr 



'n (Mt) 

t + 1 
until convergence 
output X as the matrix with enhanced form Mf. 




10 20 30 40 50 60 

Time (vectorized) 



70 80 90 100 



Figure 2. The performance of SVT for a 101 x 101 data ma- 
trix that contains 30 random frequency spikes. 5.8% of all 
entries {m = 600) are observed with signal-to-noise ampli- 
tude ratio 10. Here, n = O.lcrmax {Mt) / [4r] empirically. 
The reconstruction against the true data for the first 100 
time instances (after vectorization) are plotted. 



where n > is the soft-thresholding level^. Besides, 
in the iiT-diniensional frequency model, 'Hx°{Qt) de- 
notes the projection of Qt onto the subspace of en- 
hanced matrices (i.e. K-io\d Hankel matrices) that are 
consistent with the observed entries. Consequently, at 
each iteration, a pair (Qt, Mt) is produced by first per- 
forming singular value shrinkage and then projecting 
the outcome onto the space of K-io\d Hankel matrices 
that are consistent with observed entries. 

Fig. 2 illustrates the performance of Algorithm 1. We 
generated a 101 x 101 data matrix X with a superpo- 
sition of 30 random complex sinusoids, and revealed 
5.8% of the total entries (i.e. m = 600) uniformly at 
random. The noise was i.i.d. Gaussian-distributed giv- 
ing a signal-to-noise amplitude ratio of 10. The recon- 
structed signal is superimposed on the ground truth 
in Fig. 2. The normalized reconstruction error was 
||X-X|Jf/ |jX||p = 0.1098, validating the stability of 
EMaC in the presence of noise. 



^A good soft thresholding level is hard to pick, and 
needs to be selected by cross validation. Hence, the phase 
transisions of SVT and EMaC do not necessarily coincide. 
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5.3. Super Resolution 

The proposed EMaC algorithm works beyond tlie ran- 
dom observation model in Theorem 1. Fig. 3 con- 
siders a synthetic super resolution example, where the 
ground truth in Fig. 3 (a) contains 6 point sources with 
constant amplitude. The low-resolution observation 
in Fig. 3 (b) is obtained by measuring low-frequency 
components [— /lo, /lo] of the ground truth. Due to 
the large width of the associated point-spread func- 
tion, both the locations and amplitudes of the point 
sources are distorted in the low-resolution image. 

We apply EMaC to extrapolate high-frequency compo- 
nents up to [-/hi, /hi], where /hi//io = 2. The recon- 
struction in Fig. 3 (c) is obtained via applying directly 
inverse Fourier transform of the spectrum to avoid pa- 
rameter estimation such as the number of modes. The 
resolution is greatly enhanced from Fig. 3 (b), sug- 
gesting that EMaC is a promising approach for super 
resolution tasks. 

6. Structured Matrix Completion 

One problem closely related to our method is comple- 
tion of K-io\d Hankel matrices from a small number of 
entries. Consider for instance the 2-D model. While 
each spectrally sparse signal can be mapped to a low- 
rank 2-fold Hankel matrix, it is not clear whether all 
2-fold Hankel matrices of rank r can be written as 
the enhanced form of an object with spectral spar- 
sity r. Therefore, one can think of recovery of A'-fold 
Hankel matrices as a more general problem than the 
spectral compressed sensing problem. Indeed, Hankel 
matrix completion has found numerous applications in 
system identification (Fazel et al., 2011), natural lan- 
guage processing (Balle & Mohri, 2012), computer vi- 
sion (Sankaranarayanan et al., 2010), medical imaging 
(Shin et al., 2012), etc. 

There has been several work concerning algorithms 
and numerical experiments for Hankel matrix com- 
pletions (Fazel et al., 2003; 2011; Markovsky, 2008). 
However, to the best of our knowledge, there has 
been little theoretical guarantee that addresses di- 
rectly Hankel matrix completion. Our analysis frame- 
work in Theorem 1 can be straightforwardly adapted 
to the general K-io\d Hankel matrix completions. No- 
tice that fi2 and fi^ are defined using the SVD of X^ 
in (11) and (12), and we only need to modify the def- 
inition of /xi , as stated in the following theorem. 



Theorem 3 Consider a K-fold Hankel matrix X^ of 
rank r. The bounds in Theorem 1 and Theorem 2 con- 
tinue to hold, if the incoherence /ii is defined as the 



smallest quantity that satisfies: V (i, I) G [ni] x [712] 



max{||C7t/*A(,,,)|||,, \\A^,^i)VV*\\l] < 



nin2 



(16) 



Condition (16) requires that the left and right singular 
vectors are sufficiently uncorrelated with the observa- 
tion basis. In fact, condition (16) is a much weaker 
assumption than (10). 

It is worth mentioning that low-rank Hankel matri- 
ces can often be converted to low-rank Toeplitz coun- 
terparts. Both Hankel and Toeplitz matrices are im- 
portant forms that capture the underlying harmonic 
structures. Our results and analysis framework easily 
extend to the Toeplitz matrix completion problem. 

7. Conclusions 

We present an efficient nonparametric algorithm to es- 
timate a spectrally sparse object from its partial time- 
domain samples, which poses spectral compressed 
sensing as a low-rank Hankel structured matrix com- 
pletion problem. Under mild conditions, our algorithm 
enables recovery of the multi-dimensional unknown 
frequencies with infinite precision, which remedies the 
basis mismatch issue that arises in conventional CS 
paradigms. To the best of our knowledge, our result on 
Hankel matrix completion is also the first theoretical 
guarantee that is close to the information-theoretical 
limit (up to a logarithmic factor). 

Appendix: Proof Outline of Theorem 1 

The EMaC method has similar spirit as the well-known 
matrix completion algorithms (Candes & Recht, 2009; 
Gross, 2011), however the additional Hankel or block- 
Hankel structures on the matrices make existing the- 
oretical results inapplicable in our framework. Never- 
theless, the golfing scheme introduced in (Gross, 2011) 
lays the foundation of our analysis. We provide here a 
sketch of the proof, with detailed derivation deferred 
to (Chen & Chi, 2013) and the supplemental material. 

Denote by T the tangent space with respect to U and 
V. Let Vt be the orthogonal projection onto T, with 
Vt^ := I — Vt denoting its orthogonal complement. 
Denote by A(^ij) the orthogonal projection onto the 
subspace spanned by A/^n . The projection onto the 
space spanned by all A(-j n 's and its complement are 
defined as 



A 



E 



A 



(^,0' 



and A^=I-A. (17) 



{i,l)e[ni]x[n2] 



Suppose that CI is obtained by sampling with replace- 
ment, i.e. Q contains m indices {2^}™ j^ that are i.i.d. 
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(a) Ground truth 
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(b) Low-resolution observation (c) High-resolution reconstruction 



Figure 3. A synthetic super resolution example, where the observation (b) is taken from low-frequency components of 
the ground truth in (a), and the reconstruction (c) is done via inverse Fourier transform of extrapolated high-frequency 
components. 



generated from [ni] x [712]. We define the associated 
projection operators as An := X^i=i-4zi, and intro- 
duce another projection operator A'^ similar to An 
but with the summation only over distinct samples. 

To prove exact recovery of EMaC, it is sufficient to 
produce a dual certificate as follows. 

Lemma 2 For a location set 51 of size m, suppose that 
the sampling operator Vn obeys 



VtAVt - ^^^^^-^VrAnVT 



1 
< -. 

- 2 



// there exists a matrix W that obeys 

'{A-A'^){UV* + W)^0, 



\\Vt{W)\\ 



F - 2n^n^ ' 



(18) 



(19) 



then Xe is the unique optimizer of EMaC. 

The dual certificate is constructed via the golfing 
scheme introduced in (Gross, 2011). Specifically, we 
generate jo independent random location sets Qj ^ 
Bernoulli(g), which is sampled with replacement. Set 



q such that 1 — p = (1 — q) ° , where p 



then 



the distribution of fi is the same as fii U 112 U • • • U iljo ■ 
For small p, we have q ~ O (p/jo) through the Taylor 
expansion. Consider a small constant e < -i. The con- 
struction of a dual certificate W proceeds as follows: 

1. Set Bq = 0, and jo '■— 31ogj /g(nin2). 

2. For 1 < i < jo, let 

B, = B,_i + (9-M0, + A^) Vt {UV* - B,_i) . 

3. Set W := -{UV* -Bj„). 



We will verify that W is a valid dual certificate sat- 
isfying the conditions (19) in Lemma 2. Our con- 
struction immediately yields {A'q + A'^) {Bi) = Bi 
for 1 < J < jo- which gives 

{A-A^){UV* + W) = Q. 

Define Fi := UV* — Bi, and hence W = Fj„, one has 

Vt iF^) = {VtAVt - q'^VTAn^Vr) {F,^i) ■ 

To proceed, we present Lemma 3 which shows that An 
is sufficiently incoherent with respect to T. 

Lemma 3 There exists a constant ci > such that if 
m > CipiCsrlog{nin2) , then 
nin2. 



-VrAnVT - VtAVt 



< e 



(20) 



with probability exceeding 1 — {nin2) 

Under the assumptions of Lemma 3, we have 

\\Vt (W)||p = \\Vt (F,o)IIf < ^'"VT- < nr'i 
It remains to show \\Vt^ (^)ll !i 5 a-s follows 



Lemma 4 There exist constants C8,cg > such that 
if m > cgmayi (fiiCs, P2, fJ-a) T log {nin2), then 

VI < ^ < Jo : \\Vt^ {q'^An, + ^^) Vt {F,)\\ < 2''-^ 
with probability at least 1 — Cs{nin2)^^ ■ 
Under the assumptions of Lemma 4, we have 

jo 

\\Vt± {W)\\ < Y, W'Pt^ {q^'-^n + A^) Vt {F,)\\ < 1/2. 

So far, we have successfully shown that W is a valid 
dual certificate with high probability, and hence EMaC 
allows exact recovery with high probability. 
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